Classical simulation of quantum many-body systems with a tree tensor network 
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We show how to efficiently simulate a quantum many-body system with tree structure when its 
entanglement is bounded for any bipartite split along an edge of the tree. This is achieved by 
expanding the time-evolving block decimation simulation algorithm for time evolution from a one 
dimensional lattice to a tree graph, while replacing a matrix product state with a tree tensor network. 
As an application, we show that any one-way quantum computation on a tree graph can be efficiently 
simulated with a classical computer. 



PACS numbers: 

Introduction.— The underlying physical laws describ- 
ing classical and quantum systems are quite different, 
which makes it generally hard to efficiently simulate a 
quantum many-body evolution with a classical computer. 
However, if the quantum evolution has certain restric- 
tions, an efficient classical simulation may be possible 
[1ESSEI100- Indeed, efficient simulation algo- 
rithms are known for several restricted quantum circuit 
models 0, Q , while a general connection between entan- 
glement and the classical simulatability of quantum lat- 
tice models has been unveiled 0, 0, 0] • In one spatial di- 
mension (ID), for instance, the state of a quantum chain 
with a limited amount of entanglement between any bi- 
partition along the chain can be efficiently described us- 
ing a matrix product state (MPS) 9J. This is exploited 
by the density matrix renormalization group (DMRG) al- 
gorithm [ljj to find the ground state of the chain and by 
the time-evolving block decimation (TEBD) method |4| 
to simulate an evolution in time. Furthermore, projected 
entangled- pair states (PEPS) have been introduced as an 
extension of MPS to simulate 2D systems Q. 

In this work we consider the simulation of a quantum 
many-body system with a restricted amount of entangle- 
ment according to a tree structure. As in the ID and 2D 
approaches mentioned above, we express the d n complex 
amplitudes c^...^ of the state of n qudits (or d-level 
quantum systems), 
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in terms of a network of tensors 0, but specialize to 
the case where this tensor network (TN) has tree struc- 
ture. Given a tree TN, we explain (i) how to simulate the 
response of the system to local operations and classical 
communication (LOCC), that is to generic manipulation 
of individual qudits, including adaptive unitary transfor- 
mations and measurements and (ii) how to simulate time 
evolution. The latter is achieved by extending the TEBD 



algorithm |4|, originally proposed to simulate ID quan- 
tum lattices, so that it applies to a much broader class of 
states and physical situations. As in the ID case, the key 
to an efficient simulation is that the amount of entangle- 
ment in the system remains sufficiently bounded during 
the time evolution. 

From the perspective of theoretical computer science, 
our results imply that one-way quantum computation 
with a tree-graph cluster state can be efficiently simu- 
lated with a classical computer. One-way quantum com- 
putation with cluster states [lij . an interesting alterna- 
tive to the quantum circuit model, had previously been 
shown to be universal for quantum computation in a 2D 
lattice 14] but classically simulatable in a ID lattice fl5| . 
By extending the classical simulatability result to tree 
cluster states, we further sharpen the boundary between 
the complexities for classical and quantum computation. 

From the perspective of computational physics, our re- 
sults provide an algorithm both to find the ground state 
and to simulate time evolution in complex quantum sys- 
tems with tree structure (see also [17|), such as den- 
drimers -a class of highly branched polymers |16| . This 
algorithm, based solely on tensor multiplications and sin- 
gular value decompositions, also offers new ways to sim- 
ulate ID systems with long-range interactions. 

Canonical form of a tree TN. We use a TN with n 
open indices and tree structure, see Fig. ft|. to express 
the d n complex coefficients of the n-qudit state More 
specifically, we consider a tree network made of tensors 
with only 3 indices each |2CJ, the network therefore con- 
taining n — 2 tensors. An index connecting two tensors 
divides the network into two sub-trees A and B and the n 
qudits into two disjoint sets. We use the term bipartition 
to refer only to such divisions. The rank of an index is 
the number of values it takes. In what follows, x denotes 
the largest rank among all indices in the network. Notice 
that the tree TN depends on 0(nx 3 ) complex coefficients. 

The Schmidt decomposition of \*£) according to bipar- 




FIG. 1: (i) A tree TN for state |*> of n = 7 qudits. Each 
open index labels an orthonormal single-qudit basis, cf Eq. 
Q. Each vertex has three edges and corresponds to a tensor 
with three indices. Pairs of tensors are connected in the net- 
work through a shared index, over which there is an implicit 
summation, (ii) Canonical form of the previous tree TN. For 
each bipartition, sub-trees A and B describe Schmidt bases, 
see Eq. An empty circle on top of an edge represents a 
set of Schmidt coefficients weighting the corresponding index. 

tition A: B reads 

|¥)= Aal*^)®!*!?), ( 2 ) 

where ($k*W) = (^ ] \^) = <W, £ Q (A Q ) 2 = 1 
and where the Schmidt rank satisfies X(AB) < X- We next 
introduce the canonical form of a tree TN, which also 
consists of tensors with three indices but where each in- 
dex a shared by two tensors carries weights, see Fig. (Q. 

Definition. A tree TN is in the canonical form for bipar- 
tition A : B if (i) the weights on the connecting index a 
correspond to the Schmidt coefficients {A a } and (ii) the 
subtrees A and B describe a set of Schmidt bases {|$^)} 
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and {|$L J )}. A tree TN is in the canonical form if it is 
so for all bipartitions. 

Theorem 1. The canonical form of an n-qudit tree TN 
n can be obtained with 0(nx A ) basic operations. 
Proof: Given a bipartition A : B, |\&) can be written as 

|*)=^|^ ] )®I0L B] ), (3) 

a 

where the sets of states {J^Jl 4 ')} and {J^Jf')} for subtrees 
A and B may not be orthonormal. Our goal is to turn 
Eq. © into the Schmidt decomposition (J2J. Let 

M = XDX\ M aa , ee (<j$ l^k 41 ), (4) 

be the spectral decomposition of the matrix M of scalar 
products in A, where X'X — XX^ — I, and D TT i = 
5 TT 'd T . Then the vectors 

1^)^-1=^(1^1^]), dr>0 , (5) 

form an orthonormal set. We define X = X\fT). Anal- 
ogous considerations in B lead to an orthonormal set 
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FIG. 2: (i) Initial decomposition, Eq. (ii) Insertion of 

(X T )~ 1 X T and YY -1 , projectors onto the subspaces gen- 
erated by {|e T )} and {J/^}}- (Hi) X T Y is replaced with its 
singular value decomposition UAV, compute in 0(x 3 ) time. 
(iv) The Schmidt decomposition, Eq. is obtained after 

contracting some indices (0(\ A ) time). 

{I/17)} and a matrix Y. Eq. J2J can be rewritten as 

\^)=Y,{X T Y) Tn \e T )®\f v ). (6) 

ttj 

From the singular value decomposition (SVD) of X T Y, 

(X T Y) = UAV, (7) 

where U^U = VV^ = I and A is a diagonal matrix, we 
obtain the Schmidt coefficients A Q = A QQ and bases 

\^)=J2(U T )ar\e T }, \^ ] )=J2 V M- (8) 
r r/ 

All the above steps, summarized in Fig. take 0(\ 4 ) 
time and are repeated for each of the n — 3 bipartitions of 
the tree TN. It is not hard to see that the scalar product 
matrices M for all edges can be computed sequentially 
in an appropriate order also in time 0(n\ ) time. 

Manipulating a tree TN.— We now describe some ba- 
sic tasks involving a tree TN, that is assumed to be in 
the canonical form. First, we can compute the reduced 
density matrix for just one or two qudits. Fig. shows 
how to proceed for one qudit. 

Theorem 2. A two-qudit reduced density matrix can be 
computed with 0(md 2 x A ) basic operations, where m is 
the number of tensors in the path that connects the two 
qudits in the tree TN. [see Fig. (0}] 

Second, we can arbitrarily relocate a qudit within the 
tree TN. This is achieved by swapping the index corre- 
sponding to this qudit with other indices in the tree. 

Theorem 3. Swapping a qudit index of a tensor with 
an index of a neighboring tensor can be achieved with 
0(d 2 \ A ) basic operations. [See Fig. JJJ]. 

Third, we can update the tree TN after a unitary gate 
U has acted either on one qudit or on two neighbor- 
ing qudits. Unitarity preserves the orthogonality of the 
Schmidt bases for most bipartitions, and as a result the 
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FIG. 3: A single-qudit reduced density matrix is computed 
by contracting the TN in (i), that represents with a 

partial trace over all qudits but one. This TN reduces to the 
TN in (ii) thanks to the orthogonality of the Schmidt bases. 
In (iii) the weights have been absorbed into the three-index 
tensors (0(dx 2 ) operations) and in (iv) the remaining shared 
indices have been contracted (0(d?x 2 ) operations). 
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FIG. 4: The reduced density matrix for two qudits separated 
by m = 3 tensors is computed by contracting the TN in (i). 
This reduces to the TN in (ii) thanks to the orthogonality of 
the Schmidt bases. In (iii) the Schmidt weights have been 
absorbed into the tensors (0(dx 2 ) operations), which in turn 
are re-organized on a line made of m — 3 columns. Networks 
(iv) to (vi) illustrate how to reduce by one the number of 
columns (0(d 2 x 4 ) operations). By iteration, a single tensor 
(vii) with four open indices is obtained. 



update process involves changing only one or two ten- 
sors. When U acts on one qudit or on two qudits that 
are connected to the same three-legged tensor, that ten- 
sor simply absorbs the gate through index contraction: 

Theorem 4. Consider a two-qudit gate U acting on 
a pair of open indices of two tensors that are nearest 
neighbors in the network. The tree TN can be updated by 
replacing these two tensors, at a cost of 0(d 3 x 3 ) basic 
operations. [See Fig. JBj]. 

Efficient simulation with a tree TN.— All the above 
manipulations require computational time (and space as 
well) that scales at most linearly in the number of qudits 
n and as a small power of the maximal Schmidt rank \- 
Therefore, in those systems where the amount of entan- 
glement across all relevant bipartitions, as characterized 
by x, scales at most polynomially in n, such manipula- 
tions can be implemented efficiently. This opens up a 
number of simulation possibilities. 

For instance, we can simulate the response of the sys- 
tem to arbitrary local manipulation. Recall that LOCC 
manipulation can be decomposed as an adaptive se- 
quence of generalized local measurements mapping pure 
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FIG. 5: In order to exchange the position of indices /3 and 
7 in (i), we contract the corresponding tensors (including all 
neighboring weights) into the four-legged tensor in (ii). After 
swapping the two indices, we split the resulting tensor (iii) 
through a SVD (of a matrix of the dimension dx X dx or dx x 
X 2 ) that costs 0(d 2 x*) basic operations. The singular values 
define the new weights of the central index, whose initial rank 
X may have increased to at most dx- The old weights for 
lateral indices are detached from the two tensors to leave the 
resulting tree TN in the canonical form. 
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states into pure states. Let £ denote one such mea- 
surement on a qudit, as characterized by a set of op- 
erators {E r }, where r labels the measurement outcome. 
Outcome r occurs with probability p r = ( i S\Ef.E r \ i S'), 
in which case the state of the system becomes \^> r ) — 
E r \^)/\/Pr- To simulate £, first we randomly draw an 
outcome r according to the probability p r — tr[E r p^E}.], 
computed from the reduced density matrix of the qu- 
dit to be measured. Then a tree TN for \^f r ) is obtained 
from that of |\&) by simply absorbing operator E r into it. 
The new maximal Schmidt rank Xr satisfies Xr < X- 

An implication of the above result is that one-way 
quantum computation on a tree can be efficiently sim- 
ulated. This follows from the fact that a tree-graph clus- 
ter state has a very simple tree TN representation, with 
X = 2, while the manipulations involved in a one-way 
computation consists of LOCC. 

The simulation of a time evolution according to a two- 
body Hamiltonian H = V\ j hij is a l so possible. As in 
the case of a ID system |4| , we expand the evolution op- 
erator V — exp(-iHt) into a series of two-qudit unitary 
gates U using a Suzuki- Trotter expansion. But now, for 
each of these gates, we first bring the indices of the two 
qudits together using Theorem 3, then absorb U into the 
tree TN using Theorem 4, and finally bring the qudit 
indices back into their initial position. This generalizes 
the TEBD algorithm 0] from ID systems to a generic 
tree TN. With minimal modifications to deal with non- 
unitary gates, The TEBD algorithm can also be used 
to simulate an evolution in imaginary time according to 
V = exp(— Ht). In this way we can compute the ground 
state of H, provided H has a finite gap A > in its spec- 
trum. Recall that the expectation value of local observ- 
ables, such as the energy E = (^>\H\^) or two-point cor- 
relators, can be computed from two-qubit reduced den- 
sity matrices, which we can obtain using Theorem 2. 

As discussed in Q for ID systems, the use of a canon- 
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FIG. 6: The tensor subnetwork in (i) is contracted into the 
four-legged tensor in («), which is then decomposed into (m) 
by using a SVD (of a d\ x d\ matrix) that takes 0(d 3 x 3 ) 
basic operations. Notice that, as in Fig. ©, special attention 
is paid to first absorbing and then detaching the weights of 
the lateral indices. This guarantees that the resulting tree 
TN is in its canonical form. 



a fundamental ingredient of the simulation algorithm. 
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ical form related to the Schmidt decomposition has an 
important advantage in situations where \ is too large. 
There |'F) may still be reasonably approximated by a tree 
TN with a smaller effective \. For a given bipartition 
A : B, the best approximation to \^} is obtained by trun- 
cating decomposition @ so as to retain only the \ terms 
with greatest weight A Q . Thus, the canonical form of the 
tree TN is ideal to implement optimal bipartite trunca- 
tions. 

A tree TN may be used to simulate a many-body sys- 
tem with a genuine tree structure (determined e.g. by 
the interaction pattern) as it is the case of dendrimers 
[lij . But it can also be used to simulate a ID system 
with long-range interactions, including periodic bound- 
ary conditions. Notice that when using a MPS (that is, 
a linear TN) to represent a ID system, the number of ten- 
sors in the path connecting two random qudits scales as 
0(n) and, consequently, so does the computational time 
to simulate a gate acting on those qudits. Instead, in a 
binary tree TN, see Fig. JJJ, the path connecting any two 
qudits contains at most O(logn) tensors. This shortens 
the time required to simulate that gate by a significant 
factor log(n) jn. 




FIG. 7: (i) Tree TN for a system with tree structure, such as a 
dendrimer. (it) Two qudits in a binary tree TN are connected 
through at most 0(log(n)) tensors. (Hi) Two qudits in a linear 
TN are typically connected through 0(n) tensors. 

We conclude by noticing that in this paper we have 
explored the most general extension of the TEBD algo- 
rithm. Indeed, it appears that a tree TN — not having 
closed loops — is the most general TN to which we can 
associate a Schmidt decomposition to each of its indices, 
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A tensor with more than three indices can be decomposed (for instance, through a sequence of 

into a tree TN made of tensors with three indices each 



